One-dimensional lattice of oscillators coupled through power-law interactions: 
Continuum limit and dynamics of spatial Fourier modes 
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We study synchronization in a system of phase-only oscillators residing on the sites of a one- 
dimensional periodic lattice. The oscillators interact with a strength that decays as a power law 
of the separation along the lattice length and is normalized by a size-dependent constant. The 
exponent a of the power law is taken in the range < a < 1 . The oscillator frequency distribution 
is symmetric about its mean (taken to be zero), and is non-increasing on [0,oo). In the continuum 
limit, the local density of oscillators evolves in time following the continuity equation that expresses 
the conservation of the number of oscillators of each frequency under the dynamics. This equation 
admits as a stationary solution the unsynchronized state uniform both in phase and over the space 
of the lattice. We perform a linear stability analysis of this state to show that when it is unstable, 
different spatial Fourier modes of fluctuations have different stability thresholds beyond which they 
grow exponentially in time with rates that depend on the Fourier modes. However, numerical 
simulations show that at long times, all the non-zero Fourier modes decay in time, while only the 
zero Fourier mode (i.e., the "mean-field" mode) grows in time, thereby dominating the instability 
process and driving the system to a synchronized state. Our theoretical analysis is supported by 
extensive numerical simulations. 

PACS numbers: 05.45.Xt 
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I. INTRODUCTION: THE MODEL 

A large collection of coupled oscillating components of 
different natural frequencies getting spontaneously syn- 
chronized to oscillate at a common frequency is a phe- 
nomenon commonly observed in a variety of biological 
and physical systems, e.g., yeast cell suspensions [i|, y], 
networks of pacemaker cells in the heart 0, |j] , groups of 
flashing fireflies [a^Q , arrays of superconducting Joseph- 
son junctions [71, |8|, and many others [9[. 

The Kuramoto model provides a simple setting to in- 
vestigate the physics of synchronization [lii - ll^] . The 
model considers a large population of phase-only oscil- 
lators of distributed natural frequencies. The oscilla- 
tors are globally coupled with coupling strengths that 
are equal and independent of their spatial distribution. 
The governing dynamical equations are 



dt 
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(1) 



for i = 1, 2, . . . , A^ :^ 1. Here, 9i is a periodic variable 
of period 2tt that denotes the phase of the ith oscillator 
whose natural frequency is given by uji e R. The pa- 
rameter K > denotes the coupling strength. The fre- 
quencies are distributed according to a probability den- 
sity g(w), assumed to be one-humped and symmetric 
about its mean. It is known that above a critical cou- 
pling strength Kc = 2/(7rg(0)), the stationary state that 
the system attains at long times is a synchronized state. 
For K < Kc, however, there is no synchronization at 



long times, and each oscillator keeps running at its own 
natural frequency. 

Although the Kuramoto model has been studied exten- 
sively in recent years, much less is known about the case 
in which the oscillators are coupled with a strength that 
is a function of their spatial distribution, e.g., an inverse 
power law in the separation. Such a form of interaction 
is relevant for many situations, for example, in the study 
of long-range synchrony in a network of coupled neurons 
[l3| or in modeling flashing flreflies where the intensity 
of light signals carrying information from one firefly to 
another falls off in the three-dimensional space as the 
inverse of square of the distance from the source. 

A simple modification of the Kuramoto model ^ to 
include a power-law interaction is to consider a one- 
dimensional periodic lattice of N sites, with each site 
containing one oscillator, and with the oscillators inter- 
acting with one another with a strength that decays as a 
power law of the separation along the lattice length [ij] . 
Consequently, the governing equations are modified from 
© to 



N 
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a > 0. 



(2) 



Here, 



je is the coordinate of the jth site on the 



lattice, where e is the lattice constant. The constant N = 
X]i=i l^i ~ Xi\~"' V I is a size-dependent normalization. 
Being on a periodic lattice, we adopt the closest distance 
convention: 
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(3) 



where the total length Ne of the lattice has been taken 
to equal unity without loss of generality. From Eq. Q , 
we see that for a < 1 , the cumulative interaction of one 
oscillator with all the remaining oscillators with aligned 
phases would diverge in the continuum limit iV — > oo in 
the absence of the normalization N, which explains its 
inclusion [ij, |l5| . Note that the case a = of the model 
corresponds to the usual Kuramoto model dT]). Let us 
mention that a model with space-dependent phases as in 
Eq. ([2]), but with a different form of coupling between 
the oscillators has recently been studied in Ref. |16l |. 

A previous study of the model dH addressed by nu- 
merical simulations the issue of synchronization, in par- 
ticular, the conditions for its emergence at long times 
as a function of a for a fixed value of K [ij. It was 
reported that a critical value ttc = 2 exists, such that 
only when a < ac does the system in the limit A^ — > oo 
synchronize at a finite value of K. For a > ac, there is 
no finite value of K at which synchronization occurs. A 
recent theoretical study however predicted the value 3/2 
for ac [l4|- In Ref. [131, the model was studied without 
the normalization A^, and for finite values oi N. In this 
work, mainly numerical simulations were employed to in- 
vestigate the dependence of A^ on a and K that leads to 
synchronization. 

In this paper, we consider the model ^ with values 
of a in the range < a < 1. We take the frequency 
distribution g{uj) to be symmetric about its mean (taken 
to be zero) and non-increasing on [0, c»). We analyze 
the model ^ in the continuum limit A" — >■ (X), when the 
lattice is replaced by a continuous line characterized by 
the spatial coordinate a; G [0, 1]. Since in this limit each 
infinitesimal element dx contains a diverging number of 
oscillators, it is natural to define a local density of oscilla- 
tors p{0; uj, X, t)dx such that of all oscillators with natural 
frequency lo that are located between x and X'\-dx aX time 
t, the fraction p{6\uj,x,t)d9dx have their phase between 
and 6 + d9. Because the dynamics conserves the total 
number of oscillators with frequency uj, the time evolu- 
tion of p(6]uj,x,t) is given by the continuity equation. 
In this equation, the local densities at different spatial 
locations appear coupled due to the interaction between 
the oscillators. 

The unsynchronized state, uniform both in the phase 
6 and in the spatial coordinate x, solves the continuity 
equation in the stationary state. By performing a linear 
stability analysis of such a state, we show that when it is 
unstable, different spatial Fourier modes of fluctuations 
have different stability thresholds beyond which they 
grow exponentially in time with different rates. How- 
ever, numerical simulations show that at long times, all 
the non-zero Fourier modes decay in time, while the zero 
mode (the "mean-field" mode) grows in time and domi- 
nates the instability process, thereby driving the system 
to a synchronized state. Such a long-time dominance of 
the mean-field mode in characterizing dynamic instabil- 
ity has been observed in systems with power-law interac- 
tions evolving under a deterministic Hamiltonian dynam- 



ics [l^ . The present study demonstrates this dominance 
for the dissipative dynamics of the model ^. Our the- 
oretical predictions for the growth rates of the unstable 
modes are corroborated by numerical simulations. 

The outline of this paper is as follows. In the follow- 
ing Section, we consider the model ^ in the continuum 
limit A — > oo. and analyze the linear stability of the 
unsynchronized state. In the region in which it is un- 
stable, we derive analytic expressions for the stability 
thresholds and the growth rates of the various Fourier 
modes of fluctuations. In Section IIII[ we test our the- 
oretical predictions for the growth rates by performing 
numerical simulations, and illustrate the dominance of 
the mean-fleld mode in characterizing the long-time in- 
stability of the unsynchronized state. In simulations, we 
employ a fast numerical algorithm to compute the in- 
teraction among the oscillators, the details of which are 
relegated to the Appendix |^ In the final Section, we 
draw our conclusions. 



II. CONTINUUM LIMIT ANALYSIS 

In this section, we consider the model ^ in the con- 
tinuum limit A — > oo, and investigate in detail the issue 
of linear stability of the unsynchronized state. Finite-A 
effects are known to be quite subtle and difficult to tackle 
even for the usual Kuramoto model [ll|, [l2| , so we do not 
attempt a finite- A^ analysis of the model ^ in this work. 

In the continuum linnt, the lattice of the system 1^ 
is densely filled with sites. Let the continuous variable 
X G [0,1] stand for the spatial coordinate along the lattice 
length. Now that each infinitesimal element dx contains 
a diverging number of oscillators, we define a local den- 
sity of oscillators p{9; uj, x, t)dx such that of all oscillators 
with natural frequency oj that are located between x and 
X -\- dx &t time t, the fraction p{6; uj, x, t)d9dx have their 
phase between 9 and 9 + d9. This density is non- negative, 
27r-periodic in 9, and obeys the normalization 



2tt 



d9 p{9;uj,x,t) = 1. 



(4) 



The equations of motion ([2]) become 

d9{uj,x,t) 
dt " 

Lo + K{a)K f dffduj'dx' ""^^^^' ^\ {9'; uj' , x' , t)g{uj'), (5) 
J \x'~x\" 

where K{a) is such that N/N — > K{a) as A^ -> oo, and 
therefore, one has 



K{a)-^ 



"+1/2 dx' 



(6) 



Since we consider a in the range < a < 1, the above 
integral is convergent. The number of oscillators with fre- 
quency UJ is conserved by the dynamics, so that the time 



evolution of p[6] u, x, t) follows the continuity equation 
dp d ( d9\ .^. 

where ^ is given by Eq. ([5]). 

Next, consider the unsynchronized state uniform both 
in the phase 9 and in the spatial coordinate x: 



pa{0\uj,x,t) 



27r 



(8) 



Such a state is evidently a stationary solution of the time 
evolution ([7]). To study its linear stability under the dy- 
namics, consider adding small fluctuations, so that 

p{e;i^,x,t) = ^-+5p{9;u,x,t)] Sp{e;uj,x,t) <^1. (9) 

Ztt 

Substituting into Eq. ([7]) and keeping terms to linear 
order in 5p, we find that Sp satisfies 

dSp{9;LL!,x,t) d6p{6;uj,x,t) 

— =^ — iij- 



dt de 
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de'duj'dx' 



\x' — x\ 



5p{e'-u',x',t)g{u:'). 



(10) 



Expressing 5p in terms of its Fourier series with respect 
to the periodic variable 6 as 

oo 

5p{e-u:,x,t)^ J2 S'pk{uJ,x,t)e'''', (11) 

we find from Eq. ([T0| that only the modes fc = ±1 are 
affected by the coupling between the oscillators, and that 
Spj-i satisfies 



d6p±i{uj,x,t) 

at 



= =FitJ(5p_ti(a;,x,t) 



+!<^ /rf^'d,'^^±i(^:^,(^'). (12) 

2 7 |x'-x|" 

One thus gets a set of equations for each position a;, all 
coupled together by the second term on the right. 

Since we have a periodic lattice, it is natural to consider 
the spatial Fourier series of (5/9_|_]^(a;, cc, i): 

oo 

Kp^^{uj,x,t)= Y. Tp^,^„,{uj,t)e^^-^^. (13) 



On substituting in Eq. p2|). we get 

d5p±x^,n{^,t) 

— = T«w<5p±i^„(w,i) 



where 



c+l/2 i2TTm{x'-x) 

Am{a) = I dx' 



x-1/2 F ^1 



(15) 



Note that A™ (a) — A_m(a)- We have checked numer- 
ically that Am (a) > 0, and that it is a monotonically 
decreasing function of m, see Fi g, [l] These properties 
can also be proved analytically [201. 1^ i^ ^1^° easy to 
prove that Am(a) — >■ as to — >■ oo. Defining w = my, we 
get 

A / ^ 2 /""/^ cos(27rw) 

A,n(a) = ^TT^ / dw \ ' . (16) 

TO^ " Jq W" 



In the limit m ^- oo, the integral evaluates to a finite 
constant equal to (27r)"~^ sin(7ra/2)r(l — a), while the 
prefactor tends to zero, yielding 



lini A„j(a) = 0. 



(17) 
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FIG. 1. (Color online) Am (a) as a function of m, obtained 
by numerically evaluating the integral in Eq. (|15p . The 
different curves correspond to different a values given by 
a = 0.25, 0.5, 0.75, 0.9 (bottom to top). 



Following [ll| , we consider the right hand side of Eq. 
(ITil) to define a linear operator A as 

^'5p±i^,„(w,i) ^ 

T«W(5p±i „(a;,t) H ^ / duj (5p±i ,„(w ,t)g{uj ). 

(18) 

Then, in terms of the eigenvalues Am of the operator A, 
Eq. dHl) gives 



^P±i,m(w,i) = 5p±i,m(a;,Am)e^'"* 

with 5p^i „j satisfying 

K{a)KKjn{,a) 



(19) 



'5p±l.m(w,Am) 



2(A„i ± iijj) 



dw'(5p±im(a;',Am)5(a;'). 

(20) 

Multiplying both sides of Eq. ([20]) by .g(i^), and then 
integrating with respect to oj, one gets 



I±,m{^ — J±,m) — 0, 



(21) 



where 



I±,rn 



d^ '5p±l,m('^7^m)5(^), 



K{a)K A„i{a) 



9{^) 



duj 



(22) 
(23) 



Since /±,,„ 7^ 0, Eq. (pij) gives the dispersion relation 



J±, 



1. 



(24) 



Let us consider g(u)) to be symmetric about its mean. 
We may assume a zero mean without loss of general- 
ity (one can easily achieve this by going into a rotating 
frame of reference). Moreover, we take g{uj) to be non- 
increasing on [0,00). i.e., g{uj) > giv) whenever lo < v. 
Examples include common frequency distributions of in- 
terest like the Gaussian and the Lorentzian distributions. 
One can prove by using Theorem 2 in Ref. [21| that for 
such a g(a;), Eq. (P^ has at most one solution for A™, 
and if it exists, it is necessarily real. Then, Eq. (|24p gives 



K{a)KK„i{a) 



dw 



A. 



Ai +< 



jg{Lo) = 1. (25) 



The above equation implies that Am > 0, as otherwise the 
left hand side is negative. Since the mth spatial mode of 
fluctuations Sp^.^ „ has the time dependence Sp^i „ ^ 
gA™t ^ggg gq^ (fT9|)). it follows that this mode is either 
linearly neutrally stable corresponding to Xm = 0, or that 
it is linearly unstable corresponding to Am > 0. The 
borderline between these two behaviors is achieved at 
the critical coupling Kc , obtained by taking the limit 
Am -^ 0+ in Eq. ([25]); we get 



J^{ni) _ 



K{a)Arn{a)7rg{0)' 



(26) 



Combining the last equation with Eq. (P5|) . the growth 



rate Am > for K > Kc is given by 
2K r°° . A, 



^5(0)Ai™^ Jo >?m 



duj 



ai^) = 1- 



(27) 



With Am (a) = A_m.(a), it follows that the Fourier 
modes m and —to have the same critical thresholds and 
the same growth rates. Since A|m|(a) decreases on in- 
creasing |r7i|. we conclude that 



'(0) 



K':!'> < K''> < K, 
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< 



(28) 



The unsynchronized state ((8)) is linearly stable if all the 
spatial modes of fluctuations decay in time; using Eq. 
(PS]). this is achieved when K < Kc . On the other 
hand, it becomes unstable for K > Kc ■ 

For a = 0, when our model ([2]) reduces to the Ku- 
ramoto model ([T|), Eqs. ^ and ([TS|) give K{a)h.m{a) = 
5m,o- In this limit, the oscillators are globally cou- 
pled with equal coupling strengths, and therefore, it is 



physical to talk of only the zero Fourier mode. This 
mode has the critical coupling given by Eq. (P^ to 

be Kc ~ 2/TTg{0), such that it is linearly unstable for 
higher values of K; this matches with known results for 
the Kuramoto model [111, [i2|- For finite values of K, all 
the other modes are linearly neutrally stable because of 
diverging critical thresholds. 
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FIG. 2. (Color online) Time evolution of the observables 
ro{t),ri{t),r2{t), and r3(f) for a single realization of the ini- 
tial state {6i{0),LjJi{0);i = 1, 2, . . . , Ai"}, where the 9i's are 
chosen uniformly in [~n,n], while the uji'a are chosen from 
a Gaussian distribution with zero mean and unit variance, 
Eq. (|29p . Thus, initially, the system is in an unsynchro- 
nized state. Here, a — 0.5, and K = 15. For this value 
of a, Eqs. ^ and ^ give Ki°^ ^ 1.59577, T^i'^ « 
4.26696, kP ^ 6.53664, A"?' ^ 7.71516,..., so that the 
Fourier modes 0,1,2,3 are all linearly unstable. Conse- 
quently, ro(f), ri(t), r2(t), and r3(f) all grow in time from their 
initial values. The data in the plot are obtained from numer- 
ical simulations with TV = 2^*. 



In order to have some representative numbers for the 
critical thresholds Kc s, and for comparison with nu- 
merical simulations reported in the next section, let us 
choose a Gaussian g{Lj) with zero mean and unit vari- 
ance: 



ai^) 



1 



= --V2 



/2^ 



(29) 



Then, using Eq. (gT]), we see that for K > /d"'\ the 
growth rate Am > of the jnth mode satisfies 



K 



K 



M 



-/^Erfc 



A„ 



V2 



= 1, 



(30) 



where Erfc[a;] is the complementary error function: 
Erfc[x] ~ -y^ J^ e~* dt. In obtaining Eq. (|30p . we have 



used the result f°° dp '' ,' ./ = ^ef^^'^^Eric 
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K^J^ 
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0.05 


52.96519 


130.12769 


165.38633 


240.09713 


274.75919 


0.1 


25.89784 


59.98356 


76.10483 


107.39521 


122.97195 


0.5 


4.26696 


6.53664 


7.71516 


9.10681 


10.02918 


0.9 


1.88898 


2.04629 


2.12267 


2.18907 


2.23565 


0.95 


1.73443 


1.80472 


1.83837 


1.86673 


1.88661 



observable 



TABLE I. K, 



(m). 



1,2,. 



, 5 as a function of a for the 

(0) 



frequency distribution (|29|) . In all cases, Kc ~ 1.59577, 



III. COMPARISON WITH NUMERICAL 
SIMULATIONS 

Here, we test the theoretical predictions of the preced- 
ing Section by comparing them with results from numer- 
ical simulations. A standard procedure for simulations is 
to integrate the equations of motion ([2]). However, note 
that the equations involve computing at every integra- 
tion step a sum over N terms for each of the N oscil- 
lators, and therefore require a total computation time 
that would scale with N as N'^. In the Appendix [XI we 
show that the equations of motion can be transformed 
into a convenient form which when integrated requires 
at every integration step a computation time that scales 
with N as NlnN. We chose this alternative form of the 
equations of motion, and integrated them using a fourth- 
order Runge-Kutta algorithm with time step equal to 
0.01. The initial state of the system was chosen to have 
each oscillator phase 6 uniformly distributed in [— tt, tt] 
and frequency uj sampled from the Gaussian distribution 
(p9|) . independently of all the other oscillators. Thus, 
initially, the system is in an unsynchronized state. We 
report here simulation results for A'^ = 2^^. 

We specifically compare the theoretical and simulation 
results for the growth rates Am's. From Eq. pO|. these 

rates depend on the critical thresholds Kc'^s, and, thus, 
their comparison would indirectly provide a test of the 
theoretical predictions for Kc s given by Eqs. (|26p and 
(I29t as 



K 



-("') 



2^2 



/7rK(a)Am(a) 



(31) 



Using the above equation, we show in Table U the de- 
pendence of Kc'^s for m = 1, 2, 3, 4, 5, as a function of 
a; note that Kc is independent of a. From the table, 
we see that the thresholds for the non-zero Fourier modes 
get larger and further apart when a is close to zero. On 
the contrary, these thresholds approach the threshold for 
the zero mode when a has a value close to unity. From 
this observation, one may guess that in the limit of a ap- 
proaching zero, all thresholds excepting the one for the 
zero mode become infinitely large, while in the limit of a 
approaching unity, they all acquire the common value of 
the zero-mode threshold. 

In simulations, we monitor the time evolution of the 



N 






m = 0,1,2,, 



(32) 



Note that ro {t) corresponds to the case in which there is 
no spatial dependence of the oscillator phases, and, thus, 
characterizes the mean-field mode. 

From our analysis in the preceding Section, we see that 
for K > Kc , when the mth spatial mode of fluctuations 
is linearly unstable, rm{t) should grow exponentially in 
time: rm{t) ~ e^"**. In our simulations, we found that for 
large enough K, the observable ro(i) grows in time and 
saturates to a value of 0(1), while all the other rm(t)'s, 
after showing an initial growth in time, decay to a value 
close to zero. An illustration of such a behavior is shown 
in Fig. [5] This observation implies that in the long-time 
limit, the dynamics is dominated by the mean-field mode, 
which drives the system to a synchronized state signaled 
by ro(i) settling to a value of order 1. 

In Fig. [31 we show the temporal evolution of ro(i) and 
ri{t) for 10 different realizations of the initial state for 
three different values of a, namely, a value close to zero, a 
value close to unity, and a value in between. In each case, 
the value of K is such that the Fourier modes and 1 are 
linearly unstable; thus, K is at least larger than Kc ■ 
Consequently, ro(i) and ri(i) are expected to grow in 
time from their initial values as ri^{t) ^ e^°*,ri{t) ^ e'^^*. 

One may see from the plots for ro (t) in Fig. [31 that 
almost all realizations exhibit growth rates close to the 
theoretical value. It is then natural to average 7'o(i) over 
large enough number of initial realizations in order to 
reduce statistical fluctuations, and extract the growth 
rate Aq. We perform this exercise for a = 0.5, and display 
the result for Aq as a function of K in Fig. [H which shows 
a very good agreement with the theoretical prediction. 

From the plots of ri (i) in Fig. [31 one may note that the 
growth rates for most initial realizations are close to the 
theoretical value, although there are some realizations 
for which significant deviations are seen. On decreasing 
the value of K towards Kc , the number of this latter 
class of realizations increases, most significantly when K 



is close to K, 



(1) 



We believe that this could be due to 



pronounced finite- A^ effects close to Kc ■ Also, note from 
Eq. (|3(I|) that Ai decreases with decreasing K. Since ri{t) 
in the limit of long times decays to a close-to-zero value, 
it leaves with a very narrow time window for the growth 
ofri(t). 

From the plots in Fig. [31 one may observe that the 
growth rate of ri (i) goes to zero when a approaches zero, 
while it gets close to the rate for ro (i) when a approaches 
unity. This is in accordance with our discussion above 
on the behavior of thresholds as a approaches zero and 
unity. 



IV. CONCLUDING REMARKS 

In this paper, wc considered a system of coupled phase- 
only oscillators, in which each site of a one-dimensional 
periodic lattice of N sites contains one oscillator of dis- 
tributed natural frequency. The oscillators are coupled 
through a power-law interaction ~ {K/N)r~°' in their 
separation r along the lattice length, where the exponent 
a lies in the range < a < 1, and A^ is a size-dependent 
normalization. The oscillator frequencies are taken to be 
distributed according to a density g{uj) that is symmet- 
ric about its mean (taken to be 0) and non-increasing on 
[0,cx)). 

We studied the model in the continuum limit iV — >■ oo. 
In this limit, one may define a local density of oscilla- 
tors which evolves in time according to the continuity 
equation expressing the conservation of number of oscil- 
lators of each frequency under the dynamics. The un- 
synchronized state, uniform both in phase and in spatial 
coordinate, represents a stationary solution of the conti- 
nuity equation. By analyzing the linear stability of such 
a state, we showed that when it is unstable, different 
spatial Fourier modes of fluctuations have different sta- 
bility thresholds beyond which they grow exponentially 
in time with different rates. However, numerical simu- 
lations show that at long times, all the non-zero Fourier 
modes decay in time, while the zero mode (the "mean- 
field" mode) grows in time to dominate the instability 
process and drive the system to a synchronized state. 
Such a long-time dominance of the mean-field mode is 
known for systems with Hamiltonian systems [19|. The 
present work illustrates the dominance for dissipative dy- 
namics within the ambit of a simple model, while leaving 
the proof of its general validity as an interesting open 
issue. 
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Appendix A: A fast numerical algorithm to compute 
the interaction term in Eq. ([2]) 

Here we discuss a numerical algorithm that transforms 
the equations of motion ^ into a convenient form in 
which the interaction term is efficiently computed by a 
Fast Fourier Transform (FFT) scheme in times of order 



A^ In A'' for a system of N oscillators. Use of FFT requires 
that we choose a power of 2 for N. 

We start with Eq. ^ which we rewrite as 



^i 



K x-^ sin( 



N 






-B,) 



{d^JY 



(Al) 



where dij is the shortest distance between sites i and j 



on a periodic lattice of N sites. We set dy, 
for i 7^ J is given by 



0, while di 



\] 



if l<\]-i\< N/2, 



1 ^ ^ li ^ *l otherwise. 
Equation (jAip may be rewritten as 



(A2) 



UJ., 



K 

N 



N N 

{ cos 6i yj Vij sin Oj — sin di 2_\ ^ij 



j=i 



j=i 



(A3) 

Here, the first summation term may be interpreted as 
the zth element of the A^ x 1 column vector formed by 
the product of an A^ x A^ matrix V = \Vij\i^j=i^2,...,N 
with an A^ X 1 column vector (sin^i sin 02 • • ■ , sin^Ar)-^, 
where Vij = l/{dij)°'. Similarly, the second summation 
term may be interpreted as the ith element of the A^ x 1 
column vector formed by the product of V with an A^ x 1 
column vector (cos^i cos 6*2 . . . ,cos9]y)'^ . 

Writing out the matrix V, we see that it has the form 



V = 



Vl 
V2 



VN ... 
Vl VN 



V2 



Vl 



VN^l 
VN VN-l 



V3 


V2 




V3 




VN 


V2 


Vl_ 



(A4) 



where wi = and 

. _l i/ir~ir 



if 2 < r < N/2 + 1, 



1/(A^ - r + 1)" otherwise. 



(A5) 



Thus, F is a circulant matrix fully specified by the ele- 
ments in the first column. The remaining columns of V 
are each cyclic permutations of the elements in the first 
column, with offset equal to the column index. 
Now, note that V can be written as 



V^vil + V2P + vsP^ + ... + vnP^-\ 
where P is the N x N cyclic permutation matrix. 



(A6) 









... 1 


1 





... 







■■■ 







1 



(A7) 



Since P^ = /, the NxN identity matrix, the eigenvalues 
of P are given by Wj = e*^^(^^^)/^; j = 1,2, ... ,N, where 
Wj is the A'^-th root of unity. Equation (jA6| then imphes 
that the eigenvalues of V arc given by Aj = J2k=i '^kw'^~^ 



for j 



1,2,. 



,N. 



One may easily check that the eigenvectors of V are the 
columns of the NxN unitary discrete Fourier transform 
matrix F == -;^[fjk]j,k=i,2.,....N, where 



fjk = e-»2'rW-i)('=-i)/^ for 1 < J, k < N. (A8) 

Then, one has F^^VF ~ [AjSij]ij=i^2,...,N- 

In terms of the matrices F and F~^, one may easily 



rewrite Eq. (jA3p as 



K 



N 



N 



cose',^(P"^)yAj(Psin( 



N 



sin 6, Y^ (F- ^),j Aj (F cos ( 



(A9) 



where {Fsiiid)j (respectively, {Fcos9)j) is the 
jth clement of the A^ x 1 column vector formed 
by multiplying the matrix F with the col- 
umn vector (sin^i sin02 • • • sin^Af)-^ (respectively, 
(cos6'i cos6'2 . . . cos6'7v)"^)- Note that {Fsm9)j and 
(F cos 6)j are just discrete Fourier transforms, and may 
be computed very efficiently by standard EFT codes 
(see, e.g., Ref. [22|). The simulations reported in Section 
Hill were performed by integrating Eq. (jA9p . 
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FIG. 3. (Color online) Time evolution of ro(f) and ri(t) for 10 different realizations of the initial state {S,;(0),(jJi(0); J = 
1,2,..., A*'}, where the QiS are chosen uniformly in [— tt, tt], while the lOi^ are chosen from a Gaussian distribution with zero 
mean and unit variance, Eq. (|29[) . Thus, initially, the system is in an unsynchronized state. Each figure corresponds to a value 



Equations ((261) and ((291) 

'(2) 



give A", 

■(3) 



(1) 



25.8978, K, 



(2) 



59.9836, K, 



(3) 



of a indicated therein. For a = 0.1, we took K = 100. 

76.1048, .... For a = 0.5, we took K = 15; here, K'^'^ ^ 4.26696, K^c^' ^ 6.53664, K^'" « 7.71516, .... For a = 0.9, we took 
K = 10; here, K^^^ Ri 1.88898, ii'F' f» 2.04629, K^ ^ 2.12267, .... We have Ki°'' « 1.59577, independent of a. Thus, for these 
three values of a with the corresponding values of K, the Fourier modes and 1 are linearly unstable. As a result, in all cases, 
ro(t) and ri{t) grow in time from their initial values as ro(t) ~ e^''*,ri(t) ~ e^^'. The data in the plots are obtained from 
numerical simulations with N — 2^*. The dotted blue line in each plot shows the exponential growth with the rates Aq and Ai 
given by Eq. pO)) . 




FIG. 4. (Color online) The points denote the values of the 
growth rate Ao obtained by fitting the growth of ro (i) in time 
as ro{t) ~ e^"' for short times, while starting from an un- 
synchronized state with Gaussian g{uj) (Eq. ([29}). The data 
for ro{t) were obtained by performing numerical simulations 
with a = 0.5, A'^ = 2^^, and by averaging over 100 initial re- 
alizations. The numerical errors are of the size of the points. 
The solid line shows the theoretical curve for Ao given by Eq. 
(|30l) . From our theoretical analysis, the growth rate is zero 
for values of K below Kc , and is non-zero above. In the 
figure, the theoretical K^ is indicated. 



